I can make this simpler by doing a radial profile rather than full 2D. Of course, this assumes symmetric PSF TO do this I need:
from astropy.io import ascii, fits
import pylab as plt
%matplotlib inline
from astropy import wcs
import numpy as np
import xidplus
from xidplus import moc_routines
import pickle
#Folder containing maps
imfolder='/Users/pdh21/astrodata/COSMOS/'
pswfits=imfolder+'COSMOS_image_250_SMAP_v4.1.fits'#SPIRE 250 map
pmwfits=imfolder+'COSMOS_image_350_SMAP_v4.1.fits'#SPIRE 350 map
plwfits=imfolder+'COSMOS_image_500_SMAP_v4.1.fits'#SPIRE 500 map
#Folder containing prior input catalogue
catfolder='/Users/pdh21/astrodata/COSMOS/P3/'
#prior catalogue
prior_cat='WP5_COSMOS_XID+SPIRE_20160815.fits'
#output folder
output_folder='./'
#-----250-------------
hdulist = fits.open(pswfits)
im250phdu=hdulist[0].header
im250hdu=hdulist[1].header
im250=hdulist[1].data*1.0E3
nim250=hdulist[2].data*1.0E3
w_250 = wcs.WCS(hdulist[1].header)
pixsize250=3600.0*w_250.wcs.cd[1,1] #pixel size (in arcseconds)
hdulist.close()
#-----350-------------
hdulist = fits.open(pmwfits)
im350phdu=hdulist[0].header
im350hdu=hdulist[1].header
im350=hdulist[1].data*1.0E3
nim350=hdulist[2].data*1.0E3
w_350 = wcs.WCS(hdulist[1].header)
pixsize350=3600.0*w_350.wcs.cd[1,1] #pixel size (in arcseconds)
hdulist.close()
#-----500-------------
hdulist = fits.open(plwfits)
im500phdu=hdulist[0].header
im500hdu=hdulist[1].header
im500=hdulist[1].data*1.0E3
nim500=hdulist[2].data*1.0E3
w_500 = wcs.WCS(hdulist[1].header)
pixsize500=3600.0*w_500.wcs.cd[1,1] #pixel size (in arcseconds)
hdulist.close()
hdulist = fits.open(catfolder+prior_cat)
fcat=hdulist[1].data
hdulist.close()
inra=fcat['RA']
indec=fcat['Dec']
#---prior250--------
prior250=xidplus.prior(im250,nim250,im250phdu,im250hdu)#Initialise with map, uncertianty map, wcs info and primary header
prior250.prior_cat(inra,indec,prior_cat)#Set input catalogue
prior250.prior_bkg(-5.0,5)#Set prior on background (assumes Guassian pdf with mu and sigma)
#---prior350--------
prior350=xidplus.prior(im350,nim350,im350phdu,im350hdu)
prior350.prior_cat(inra,indec,prior_cat)
prior350.prior_bkg(-5.0,5)
#---prior500--------
prior500=xidplus.prior(im500,nim500,im500phdu,im500hdu)
prior500.prior_cat(inra,indec,prior_cat)
prior500.prior_bkg(-5.0,5)
#pixsize array (size of pixels in arcseconds)
pixsize=np.array([pixsize250,pixsize350,pixsize500])
#point response function for the three bands
prfsize=np.array([18.15,25.15,36.3])
#use Gaussian2DKernel to create prf (requires stddev rather than fwhm hence pfwhm/2.355)
from astropy.convolution import Gaussian2DKernel
##---------fit using Gaussian beam-----------------------
prf250=Gaussian2DKernel(prfsize[0]/2.355,x_size=101,y_size=101)
prf250.normalize(mode='peak')
prf350=Gaussian2DKernel(prfsize[1]/2.355,x_size=101,y_size=101)
prf350.normalize(mode='peak')
prf500=Gaussian2DKernel(prfsize[2]/2.355,x_size=101,y_size=101)
prf500.normalize(mode='peak')
pind250=np.arange(0,101,1)*1.0/pixsize[0] #get 250 scale in terms of pixel scale of map
pind350=np.arange(0,101,1)*1.0/pixsize[1] #get 350 scale in terms of pixel scale of map
pind500=np.arange(0,101,1)*1.0/pixsize[2] #get 500 scale in terms of pixel scale of map
prior250.set_prf(prf250.array,pind250,pind250)#requires psf as 2d grid, and x and y bins for grid (in pixel scale)
prior350.set_prf(prf350.array,pind350,pind350)
prior500.set_prf(prf500.array,pind500,pind500)
order=10
Tile=6978147
moc=moc_routines.get_fitting_region(order,Tile)
prior250.set_tile(moc)
prior350.set_tile(moc)
prior500.set_tile(moc)
prior250.get_pointing_matrix_unknown_psf()
prior350.get_pointing_matrix()
prior500.get_pointing_matrix()
import pylab as plt
%matplotlib inline
prior250.upper_lim_map()
prior350.upper_lim_map()
prior500.upper_lim_map()
print prior250.nsrc
from xidplus.stan_fit import SPIRE
fit=SPIRE.all_bands_PSF(prior250,prior350,prior500,51,np.arange(51),iter=1000)
posterior=xidplus.posterior_stan(fit,[prior250,prior350,prior500])
for i in range(0,2000):
plt.plot(posterior.stan_fit[:,:,-52:-1].reshape(2000,51)[i,:],'b',alpha=0.02,linestyle='-')
plt.plot()
plt.plot(prior250.prf[51:,51],'r')
plt.figure(figsize=(10,5))
psf_fit=np.percentile(posterior.stan_fit[:,:,-52:-1].reshape(2000,51),[16.0,50.0,84.0],axis=0)
plt.plot(psf_fit[1,:],'b')
plt.fill_between(np.arange(0,51),psf_fit[0,:],psf_fit[2,:],alpha=0.2)
plt.plot(prior250.prf[51:,51],'r')
plt.xlabel(r'$R (arcsec.)$')
plt.ylabel(r'PSF normalisation')
Sigma=np.empty((51,51))
for i in range(0,51):
for j in range(0,51):
Sigma[i,j] =np.exp(-0.01 * np.power(np.arange(51)[i] - np.arange(51)[j],2))
if i==j:
Sigma[i,j]=Sigma[i,j]+0.00001
mu=np.zeros(51)
for i in range(0,10):
plt.plot(np.random.multivariate_normal(mu, Sigma,10)[i])
def ymod_map(prior,flux,PSF):
from scipy.sparse import coo_matrix
f=coo_matrix((flux, (range(0,prior.nsrc),np.zeros(prior.nsrc))), shape=(prior.nsrc, 1))
A=coo_matrix((PSF[prior.amat_data.astype(long)], (prior.amat_row, prior.amat_col)), shape=(prior.snpix, prior.nsrc))
rmap_temp=(A*f)
return np.asarray(rmap_temp.todense())
from xidplus import posterior_maps as postmaps
hdulist_250=postmaps.make_fits_image(prior250,prior250.sim)
hdulist_350=postmaps.make_fits_image(prior350,prior350.sim)
hdulist_500=postmaps.make_fits_image(prior500,prior500.sim)
#Set some color info
import seaborn as sns
cmap=sns.cubehelix_palette(8, start=.5, rot=-.75,as_cmap=True)
vmin=-1.7E1
vmax=4.446e+01
import warnings
from matplotlib.cbook import MatplotlibDeprecationWarning
warnings.simplefilter('ignore', MatplotlibDeprecationWarning)
warnings.simplefilter('ignore', UserWarning)
warnings.simplefilter('ignore', RuntimeWarning)
warnings.simplefilter('ignore',UnicodeWarning)
import aplpy
fig = plt.figure(figsize=(30,10))
real_250 = aplpy.FITSFigure(hdulist_250[1],figure=fig,subplot=(1,3,1))
real_250.show_colorscale(vmin=vmin,vmax=vmax,cmap=cmap)
real_250.show_markers(prior250.sra, prior250.sdec, edgecolor='black', facecolor='black',
marker='o', s=20, alpha=0.5)
real_250.tick_labels.set_xformat('dd.dd')
real_250.tick_labels.set_yformat('dd.dd')
real_350 = aplpy.FITSFigure(hdulist_350[1],figure=fig,subplot=(1,3,2))
real_350.show_colorscale(vmin=vmin,vmax=vmax,cmap=cmap)
real_350.show_markers(prior250.sra, prior250.sdec, edgecolor='black', facecolor='black',
marker='o', s=20, alpha=0.5)
real_350.tick_labels.set_xformat('dd.dd')
real_350.tick_labels.set_yformat('dd.dd')
real_500 = aplpy.FITSFigure(hdulist_500[1],figure=fig,subplot=(1,3,3))
real_500.show_colorscale(vmin=vmin,vmax=vmax,cmap=cmap)
real_500.show_markers(prior250.sra, prior250.sdec, edgecolor='black', facecolor='black',
marker='o', s=20, alpha=0.5)
real_500.tick_labels.set_xformat('dd.dd')
real_500.tick_labels.set_yformat('dd.dd')
mod_map=np.full((hdulist_250[1].data.shape[1],hdulist_250[1].data.shape[0],500),np.nan)
mod_map_array=np.empty((prior250.snpix,500))
for i in range(0,500):
mod_map_array[:,i]= ymod_map(prior250,posterior.stan_fit[i,0,0:prior250.nsrc],posterior.stan_fit[i,0,-52:-1]).reshape(-1)+posterior.stan_fit[i,0,prior250.nsrc]+np.random.normal(scale=np.sqrt(prior250.snim**2+posterior.stan_fit[i,0,(prior250.nsrc+1)*3]**2))
mod_map[prior250.sx_pix-np.min(prior250.sx_pix)-1,prior250.sy_pix-np.min(prior250.sy_pix)-1,i]=mod_map_array[:,i]
mod_map_350=np.full((hdulist_350[1].data.shape[1],hdulist_350[1].data.shape[0],500),np.nan)
mod_map_array_350=np.empty((prior350.snpix,500))
for i in range(0,500):
mod_map_array_350[:,i]= postmaps.ymod_map(prior350,posterior.stan_fit[i,0,prior350.nsrc+1:2*prior350.nsrc+1]).reshape(-1)+posterior.stan_fit[i,0,2*prior350.nsrc+1]+np.random.normal(scale=np.sqrt(prior350.snim**2+posterior.stan_fit[i,0,1+(prior350.nsrc+1)*3]**2))
mod_map_350[prior350.sx_pix-np.min(prior350.sx_pix)-1,prior350.sy_pix-np.min(prior350.sy_pix)-1,i]=mod_map_array_350[:,i]
mod_map_500=np.full((hdulist_500[1].data.shape[1],hdulist_500[1].data.shape[0],500),np.nan)
mod_map_array_500=np.empty((prior500.snpix,500))
for i in range(0,500):
mod_map_array_500[:,i]= postmaps.ymod_map(prior500,posterior.stan_fit[i,0,2*prior500.nsrc+2:3*prior350.nsrc+2]).reshape(-1)+posterior.stan_fit[i,0,3*prior500.nsrc+2]+np.random.normal(scale=np.sqrt(prior500.snim**2+posterior.stan_fit[i,0,2+(prior500.nsrc+1)*3]**2))
mod_map_500[prior500.sx_pix-np.min(prior500.sx_pix)-1,prior500.sy_pix-np.min(prior500.sy_pix)-1,i]=mod_map_array_500[:,i]
import copy
prior250_old=copy.copy(prior250)
prior250_old.get_pointing_matrix()
plt.plot(mod_map_array[:,0]/mod_map_array_2[:,0])
from tempfile import NamedTemporaryFile
VIDEO_TAG = """<video controls>
<source src="data:video/x-m4v;base64,{0}" type="video/mp4">
Your browser does not support the video tag.
</video>"""
def anim_to_html(anim):
if not hasattr(anim, '_encoded_video'):
with NamedTemporaryFile(suffix='.mp4') as f:
anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])
video = open(f.name, "rb").read()
anim._encoded_video = video.encode("base64")
return VIDEO_TAG.format(anim._encoded_video)
from matplotlib import animation
from IPython.display import HTML
def display_animation(anim):
plt.close(anim._fig)
return HTML(anim_to_html(anim))
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=(30,10))
res250=aplpy.FITSFigure(hdulist_250[1],figure=fig,subplot=(1,3,1))
res250.show_colorscale(vmin=vmin,vmax=vmax,cmap=cmap)
res250.show_markers(prior250.sra, prior250.sdec, edgecolor='black', facecolor='black',
marker='o', s=20, alpha=0.5)
res350=aplpy.FITSFigure(hdulist_350[1],figure=fig,subplot=(1,3,2))
res350.show_colorscale(vmin=vmin,vmax=vmax,cmap=cmap)
res350.show_markers(prior350.sra, prior350.sdec, edgecolor='black', facecolor='black',
marker='o', s=20, alpha=0.5)
res500=aplpy.FITSFigure(hdulist_500[1],figure=fig,subplot=(1,3,3))
res500.show_colorscale(vmin=vmin,vmax=vmax,cmap=cmap)
res500.show_markers(prior500.sra, prior500.sdec, edgecolor='black', facecolor='black',
marker='o', s=20, alpha=0.5)
# animation function. This is called sequentially
def animate(i):
res250._data=mod_map[:,:,i].T
res350._data=mod_map_350[:,:,i].T
res500._data=mod_map_500[:,:,i].T
res250.show_colorscale(vmin=vmin,vmax=vmax,cmap=cmap)
res350.show_colorscale(vmin=vmin,vmax=vmax,cmap=cmap)
res500.show_colorscale(vmin=vmin,vmax=vmax,cmap=cmap)
return [res250,res350,res500]
# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate,
frames=500, interval=1000)
anim.save('post_animation.mp4', fps=20, extra_args=['-vcodec', 'libx264'])
# call our new function to display the animation
display_animation(anim)